********************************************************************************
********************************************************************************
clear
set more off
* change these paths depending on where the dataset is saved and where you
* want to generate the output
*global datapath "~/Dropbox/EmpireGIS/Paper/Replication"
*global resultspath "~/Dropbox/EmpireGIS/Paper/Replication"


********************************************************************************
*  File-Name:       turnovers_analysis_Replication.do
*  Date:            20 Jan 2019
*  Author:          S.P.Harish & Christopher Paik
*  Purpose:         Replication for the Historical State Stability and Economic 
*                   Development in Europe
*  Previous File:   None
*  Data Used:       turnovers_replication_20191020.dta
*  Output File:     Tables & Figures as in the paper
*  O/S:             Mac
*  Software Needed: esttab (for table output)
*                   Note: you need to install the spatial program to run the 
*                   below code
*                   ssc install reg2hdfe
*                   ssc install tmpdir
*                   
*                   Note: you will need to manually copy the ado files for 
*                   ols_spatial_HAC and reg2hdfespatial to your ado folder
********************************************************************************




********************************************************************************
* Start with the master file
********************************************************************************

use "$datapath/turnovers_replication_20191021.dta", clear 


********************************************************************************
* Table 1 
********************************************************************************

eststo clear

estpost summarize gdp2010percapita_maxarea gdp2010percapita_avg gdp2010percapita_weighted ///
 iv2000a SQiv2000a ///
 altiv2000a SQaltiv2000a ///
altivTV2000a SQaltivTV2000a ///
anystate2000a SQanystate2000a ///
agricdate ramankutty_mean riverdist elevation_mean dist2city ///
romanempire ottomanempire mongolianempire lat lon latlon ///
east_europe low_region westgermany eubelt capital num_centuries ///
mean_parl mean_plunder ///
num_siege num_battle ///
if flag ==1 



esttab using "$resultspath/SummaryStats.tex"	///
	, ///
	cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") ///
	substitute(m$ m\\$)	///
	nonumber label replace unstack nomtitle noobs

eststo clear



********************************************************************************
* Table 2
********************************************************************************

foreach dv of varlist loggdp2010percapita_maxarea {

eststo clear

eststo: xi: reg `dv' iv2000a SQiv2000a i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean elevation_mean i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean elevation_mean riverdist i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean elevation_mean riverdist agricdate i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean elevation_mean riverdist agricdate dist2city i.entity_id2000, robust
eststo: xi: reg `dv' iv2000a SQiv2000a ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon i.entity_id2000, robust

esttab using "$resultspath/`dv'.tex", booktabs label replace ///
compress nomtitles nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
indicate("Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = lat" "Longitude = lon" "Lat*Lon = latlon" "Country FE in Yr2000 = _I*", label($\checkmark$ ))	///
drop(_cons) nonotes

eststo clear

}



********************************************************************************
* Figures 5 and 6
********************************************************************************

xi: reg loggdp2010percapita_maxarea iv2000a c.iv2000a##c.iv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

margins, at(iv2000a==(0(0.1)1))
marginsplot, title("Predicted Effect of Mean State Duration") ///
ytitle("2010 GDP per capita (Log)", size (medsmall) )  ///
xtitle("Mean State Duration", size(medsmall))  ///
ylabel(, labsize(small))  ///
xlabel(, labsize(small))  ///
yline(0)  ///
scheme(s1mono)

graph export "$resultspath/QEmsd.png", as(png) replace

xi: reg loggdp2010percapita_maxarea iv2000a c.iv2000a##c.iv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

margins, dydx(iv2000a) at(iv2000a==(0(0.1)1))
marginsplot, title("Marginal Effect of Mean State Duration") ///
ytitle("Change in 2010 GDP per capita (Log)", size (medsmall) )  ///
xtitle("Mean State Duration", size(medsmall))  ///
ylabel(, labsize(small))  ///
xlabel(, labsize(small))  ///
yline(0)  ///
scheme(s1mono)

graph export "$resultspath/MEmsd.png", as(png) replace





********************************************************************************
* Table 3
********************************************************************************

eststo clear

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
anystate2000a SQanystate2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
romanempire ottomanempire mongolianempire ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
east_europe low_region westgermany eubelt capital ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
east_europe low_region westgermany eubelt num_centuries ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust


esttab using "$resultspath/robust1.tex", booktabs label replace ///
compress nomtitles nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
indicate("State Presence = anystate2000a" "State Presence Sq = SQanystate2000a" ///
"Roman Empire = romanempire" "Ottoman Empire = ottomanempire" "Mongolian Empire = mongolianempire" ///
"Eastern Europe = east_europe" "Low Region = low_region" "West Germany = westgermany" ///
"European City Belt = eubelt" "European Capital Cities = capital" ///
 "Number of Centuries with Capital = num_centuries" ///
"Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = lat" "Longitude = lon" "Lat*Lon = latlon" "Country FE in Yr2000 = _I*", label($\checkmark$ ))	///
drop(_cons) nonotes

eststo clear



********************************************************************************
* Table 4
********************************************************************************

eststo clear

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
mean_parl  ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
mean_plunder ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
mean_parl mean_plunder ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
num_siege  ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
num_battle  ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_maxarea iv2000a SQiv2000a ///
num_siege num_battle ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust


esttab using "$resultspath/robust2.tex", booktabs label replace ///
compress nomtitles nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
indicate("Parliament = mean_parl" "Plunder = mean_plunder" ///
"Number of Sieges = num_siege" "Number of Battles = num_battle" ///
"Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = lat" "Longitude = lon" "Lat*Lon = latlon" "Country FE in Yr2000 = _I*", label($\checkmark$ ))	///
drop(_cons) nonotes

eststo clear



********************************************************************************
* Table 5
********************************************************************************

eststo clear

eststo: xi: reg  loggdp2010percapita_maxarea altiv2000a SQaltiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg  loggdp2010percapita_maxarea  altivTV2000a SQaltivTV2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_avg iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

eststo: xi: reg loggdp2010percapita_weighted iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

esttab using "$resultspath/robust3.tex", booktabs label replace ///
compress nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("\shortstack{GDPPC\\Max Area}" "\shortstack{GDPPC\\Max Area}" "\shortstack{GDPPC\\Average}" "\shortstack{GDPPC\\Weighted}") ///
indicate("Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = lat" "Longitude = lon" "Lat*Lon = latlon" "Country FE in Yr2000 = _I*", label($\checkmark$ ))	///
drop(_cons) nonotes

eststo clear





*******************************************************************************
* Figure S1
*******************************************************************************

local years1 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100
foreach i of local years1 {

xi: reg loggdp2010percapita_maxarea msd`i'b c.msd`i'b##c.msd`i'b ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

margins, at(msd`i'b==(0(0.1)1)) 
marginsplot, title("`i'-2000", size(small)) ///
ytitle("")  ///
xtitle("")  ///
scheme(s1mono)

graph save "$resultspath/QEmsd`i'b.gph", replace
*graph export "../TablesFigures/QEmsd`i'b.png", as(png) replace

}

graph combine ///
"$resultspath/QEmsd1900b.gph" ///
"$resultspath/QEmsd1800b.gph" ///
"$resultspath/QEmsd1700b.gph" ///
"$resultspath/QEmsd1600b.gph" ///
"$resultspath/QEmsd1500b.gph" ///
"$resultspath/QEmsd1400b.gph" ///
"$resultspath/QEmsd1300b.gph" ///
"$resultspath/QEmsd1200b.gph" ///
"$resultspath/QEmsd1100b.gph" ///
"$resultspath/QEmsd1000b.gph" ///
"$resultspath/QEmsd900b.gph" ///
"$resultspath/QEmsd800b.gph" ///
"$resultspath/QEmsd700b.gph" ///
"$resultspath/QEmsd600b.gph" ///
"$resultspath/QEmsd500b.gph" ///
"$resultspath/QEmsd400b.gph" ///
"$resultspath/QEmsd300b.gph" ///
"$resultspath/QEmsd200b.gph" ///
"$resultspath/QEmsd100b.gph" ///
, ///
title("Quadratic Relationship Over Time") ///
l1title("Log 2010 GDP per capita", size (medsmall) )  ///
b1title("Mean State Duration", size(medsmall))  ///
rows(3) cols (7) scheme(s1mono)

graph export "$resultspath/QEmsdComb.png", as(png) replace

***********************
* Figure 7
***********************

graph combine ///
"$resultspath/QEmsd500b.gph" ///
"$resultspath/QEmsd800b.gph" ///
"$resultspath/QEmsd1400b.gph" ///
, ///
title("Quadratic Relationship Over Time") ///
l1title("Log 2010 GDP per capita", size (medsmall) )  ///
b1title("Mean State Duration", size(medsmall))  ///
rows(3) cols (7) scheme(s1mono)

graph export "$resultspath/QEmsd3period.png", as(png) replace



*******************************************************************************
* Figure S2
*******************************************************************************


local years1 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100
foreach i of local years1 {

xi: reg loggdp2010percapita_maxarea msd`i'a c.msd`i'a##c.msd`i'a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

margins, at(msd`i'a==(0(0.1)1)) 
marginsplot, title("0000-`i'", size(small)) ///
ytitle("")  ///
xtitle("")  ///
scheme(s1mono)

graph save "$resultspath/QEmsd`i'a.gph", replace
*graph export "../TablesFigures/QEmsd`i'b.png", as(png) replace

}

graph combine ///
"$resultspath/QEmsd1900a.gph" ///
"$resultspath/QEmsd1800a.gph" ///
"$resultspath/QEmsd1700a.gph" ///
"$resultspath/QEmsd1600a.gph" ///
"$resultspath/QEmsd1500a.gph" ///
"$resultspath/QEmsd1400a.gph" ///
"$resultspath/QEmsd1300a.gph" ///
"$resultspath/QEmsd1200a.gph" ///
"$resultspath/QEmsd1100a.gph" ///
"$resultspath/QEmsd1000a.gph" ///
"$resultspath/QEmsd900a.gph" ///
"$resultspath/QEmsd800a.gph" ///
"$resultspath/QEmsd700a.gph" ///
"$resultspath/QEmsd600a.gph" ///
"$resultspath/QEmsd500a.gph" ///
"$resultspath/QEmsd400a.gph" ///
"$resultspath/QEmsd300a.gph" ///
"$resultspath/QEmsd200a.gph" ///
"$resultspath/QEmsd100a.gph" ///
, ///
title("Quadratic Relationship Over Time") ///
l1title("Log 2010 GDP per capita", size (medsmall) )  ///
b1title("Mean State Duration", size(medsmall))  ///
rows(3) cols (7) scheme(s1mono)

graph export "$resultspath/QEmsdCombA.png", as(png) replace


local years1 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100
foreach i of local years1 {

rm "$resultspath/QEmsd`i'a.gph" 
rm "$resultspath/QEmsd`i'b.gph"  

}






*******************************************************************************
* Table S1
*******************************************************************************

eststo clear

foreach var of varlist ppltrst revipeqopt revimpfree ipfrule {

eststo: xi: reg `var' iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon i.entity_id2000, robust

}


esttab using "$resultspath/ess.tex", booktabs label replace ///
compress mtitles("Trust" "Respect" "Control" "Obedience") nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
indicate("Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = lat" "Longitude = lon" "Lat*Lon = latlon" "Country FE in Yr2000 = _I*", label($\checkmark$ ))	///
drop(_cons) nonotes

eststo clear


*******************************************************************************
* Figure S3
*******************************************************************************

local years1 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100
foreach i of local years1 {

xi: reg revimpfree msd`i'a c.msd`i'a##c.msd`i'a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city lat lon latlon ///
i.entity_id2000, robust

margins, at(msd`i'a==(0(0.1)1)) 
marginsplot, title("0000-`i'", size(small)) ///
ytitle("")  ///
xtitle("")  ///
scheme(s1mono)

graph save "$resultspath/QEctrl`i'a.gph", replace
*graph export "../TablesFigures/QEmsd`i'b.png", as(png) replace

}

graph combine ///
"$resultspath/QEctrl1900a.gph" ///
"$resultspath/QEctrl1800a.gph" ///
"$resultspath/QEctrl1700a.gph" ///
"$resultspath/QEctrl1600a.gph" ///
"$resultspath/QEctrl1500a.gph" ///
"$resultspath/QEctrl1400a.gph" ///
"$resultspath/QEctrl1300a.gph" ///
"$resultspath/QEctrl1200a.gph" ///
"$resultspath/QEctrl1100a.gph" ///
"$resultspath/QEctrl1000a.gph" ///
"$resultspath/QEctrl900a.gph" ///
"$resultspath/QEctrl800a.gph" ///
"$resultspath/QEctrl700a.gph" ///
"$resultspath/QEctrl600a.gph" ///
"$resultspath/QEctrl500a.gph" ///
"$resultspath/QEctrl400a.gph" ///
"$resultspath/QEctrl300a.gph" ///
"$resultspath/QEctrl200a.gph" ///
"$resultspath/QEctrl100a.gph" ///
, ///
title("Quadratic Relationship Over Time") ///
l1title("2010 Control Levels", size (medsmall) )  ///
b1title("Mean State Duration", size(medsmall))  ///
rows(3) cols (7) scheme(s1mono)

graph export "$resultspath/QEctrlCombA.png", as(png) replace


local years1 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 700 600 500 400 300 200 100
foreach i of local years1 {

rm "$resultspath/QEctrl`i'a.gph"

}



********************************************************************************
* Table 6
********************************************************************************

* generate constant time var
cap drop year 
cap drop _merge
gen year = 2000

* rename dv bcoz of length
ren loggdp2010percapita_maxarea dv
* rename lat lon bcoz they are reserved in conley code
ren lat y
ren lon x
* create copies of lat and lon bcoz cannot use the same var in conley code
gen latitude = y
gen longitude = x

eststo clear

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
, ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean , ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean elevation_mean , ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist , ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate , ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city , ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)

eststo: reg2hdfespatial dv iv2000a SQiv2000a ///
ramankutty_mean elevation_mean riverdist agricdate dist2city latlon latitude longitude, ///
timevar(year) panelvar(entity_id2000) lat(y) lon(x) distcutoff(200) lagcutoff(0)


esttab using "$resultspath/robust4.tex", booktabs label replace ///
compress nomtitles nogaps ///
nodepvars se(3) b(3) ///
stats(N, labels(Observations) fmt(0 3))  ///
star(* 0.10 ** 0.05 *** 0.01) ///
indicate("Agricultural Suitability = ramankutty_mean" "Elevation = elevation_mean" ///
"Distance to Water = riverdist" "Agri Adoption = agricdate" "Distance to City = dist2city" ///
"Latitude = latitude" "Longitude = longitude" "Lat*Lon = latlon" , label($\checkmark$ ))	///
nonotes


eststo clear



STOP
